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ABSTRACT 


Under  an  age  replacement  policy  a  system  is  replaced  at  a  fixed  age  q>  or  at  failure 
whichever  comes  first.  If  the  cost  of  replacing  the  system  before  failure  is  less  than  the 
cost  of  replacing  it  at  failure,  this  type  of  maintenance  policy  can  lead  to  considerable 
savings.  An  often  used  criterion  for  finding  an  "optimal"  replacement  age  <p,  is  to  mini¬ 
mize  the  long  run  expected  cost  per  unit  time  of  a  policy  with  replacement  age  cp.  This 
cost  function  clearly  depends  on  the  underlying  distribution  of  the  system  lifetimes. 
When  this  distribution  is  unknown,  the  cost  function  and  hence  <p*  need  to  be  estimated. 

In  this  thesis,  we  study  the  large  and  small  sample  properties  of  a  procedure  which 
estimates  <p*.  In  particular,  we  study  sequential  maximum  likelihood  estimators  of  <p* 
which  are  updated  at  each  replacement  based  on  the  replacement  history  of  the  system 
so  far.  In  this  sequential  procedure  each  system  is  subject  to  the  age  replacement  policy 
with  estimated  <p*  based  on  all  the  data  gathered  so  far.  This  type  of  procedure  should 
control  the  actual  cost  per  unit  time  while  gathering  data  needed  to  estimate  q>*. 

This  thesis  contains  a  detailed  description  of  the  sequential  estimation  procedure 
when  the  underhing  system  life  times  have  a  Weibull  distribution  and  a  Gamma  dis¬ 
tribution.  Monte-Carlo  methods  are  then  used  to  study  the  behavior  of  the  estimated 
optimal  age  replacement  times  and  more  importantly  the  actual  costs  per  unit  time  for 
different  sample  sizes,  costs  and  choices  of  the  underhing  Weibull  and  Gamma  distrib¬ 
utions. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Optimal  maintenance  policies  are  designed  to  reduce  the  number  of  system  failures 
and  minimize  the  cost  of  repair  by  scheduling  planned  replacements.  By  far,  most  of  the 
research  in  this  area  has  been  from  the  modeling  stand-point.  Even  in  the  most  basic 
scenario,  where  the  underlying  system  lifetimes  are  assumed  to  be  independent  and 
identically  distributed,  the  problem  of  updating  the  maintenance  policy  using  the  past 
maintenance  history  has  not  been  adequately  solved.  In  this  thesis,  Monte-Carlo 
methods  are  used  to  study  a  particular  parametric  procedure  which  updates  estimates 
of  an  optimal  age  replacement  policy  after  each  replacement. 

B.  MAINTENANCE  AND  MAINTENANCE  POLICIES 

Maintenance  is  a  combination  of  actions  carried  out  to  retain  a  unit  in,  or  restore 
it  to.  an  acceptable  condition.  There  are  two  forms  of  maintenance.  These  are  preventive 
and  corrective  maintenance.  Corrective  maintenance  is  carried  cut  to  restore  (including 
minor  adjustment  and  repair)  a  unit  which  has  ceased  to  meet  an  acceptable  condition. 
Preventive  maintenance  is  interned  to  reduce  the  likelihood  of  a  unit  not  meeting  an 
acceptable  condition  [Ref.  1:  pp.  1-8]. 

Maintenance  policies  reduce  the  number  of  system  failures  and  maintenance  costs 
by  adopting  a  schedule  of  planned  maintenance.  For  instance,  when  the  failure  of  a  unit 
during  actual  operation  is  dangerous  or  costly  (i.e.,  the  cost  C,  of  unscheduled  mainte¬ 
nance  due  to  system  failure  is  more  than  the  cost  C2  of  scheduled  maintenance  before 
system  failure)  and  the  unit  is  characterized  by  a  failure  rate  that  increases  with  age.  it 
may  be  wise  to  maintain  it  before  it  has  aged  too  greatly  (Ref.  2:  p.  46].  On  the  other 
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hand,  performing  maintenance  actions  too  frequently  can  also  be  costly.  Thus,  choosing 
the  preventive  maintenance  policy  that  schedules  maintenance  to  best  control  costs  will 
be  influenced  by  the  relative  costs  of  failure  and  preventive  maintenance  and  the 
stochastic  properties  of  the  lifetime  of  the  unit  [Ref.  3:  pp.  19-24].  We  are  interested 
in  determining  the  sequence  of  times  at  which  preventive  maintenance  should  occur.  In 
particular,  we  are  actually  interested  in  determining  an  optimal  maintenance  policies 
which  minimizes  long  run  expected  maintenance  cost  per  unit  time. 

Preventive  maintenance  is  the  total  of  all  service  functions  aimed  at  maintaining  and 
improving  reliability  performance  characteristics  and  concerns  itself  with  such  activ  ities 
as  the  replacement  and  repair  of  systems,  inspections,  testing  and  checking  of  working 
parts  during  their  operation.  In  this  thesis,  we  will  only  consider  maintenance  actions 
which  involve  replacement  of  a  system.  It  will  be  assumed  that  the  replacement  action 
returns  the  equipment  to  the  as  new  condition,  thus  providing  the  same  services  as  the 
equipment  which  has  just  been  replaced.  By  making  this  assumption,  we  arc  implying 
that  the  distribution  of  time  to  failure  of  the  new  svstem  is  the  same  as  that  of  the  system 
which  was  replaced.  In  addition,  we  assume  that  system  lifetimes  are  independent.  The 
unscheduled  and  scheduled  replacement  costs  C,  and  C2  remain  constant  where  C,  >  C2. 

An  important  replacement  policy  is  the  policy  based  on  age  (age  replacement).  Such 
a  police  is  in  force  if  a  unit  is  always  replaced  at  the  time  of  failure  or  (p  units  of  time 
after  its  installation,  whichever  comes  first.  Under  a  block  replacement  policy  the  unit 
is  replaced  at  times  k<p  (k  =  1,2, ...),  and  at  failure.  This  replacement  policy  derives  its 
name  from  the  commonly  employed  practice  of  replacing  a  block  or  group  of  units  in  a 
system  at  prescribed  times  k<p  (k  =  1,2,  ...)  independent  of  the  failure  h’Story  of  the 
system  [Ref.  2:  p.  46]. 
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"  Age  replacement  is  administratively  more  difficult  to  implement,  since  the  age  of 
the  unit  must  be  recorded.  But  block  replacement,  although  simpler  to  administer  since 
the  age  of  the  system  need  not  be  recorded,  leads  to  more  frequent  replacement  of  rela¬ 
tively  new  items  "  [Ref.  4:  p.  158].  In  this  thesis  we  only  consider  age  replacement 
policies. 

An  optimal  age  replacement  policy  is  the  age  replacement  policy  which  yields  the 
smallest  long  run  expected  replacement  costs  per  unit  time.  To  find  an  optimal  age  re¬ 
placement  policy  we  require  explicit  knowledge  about  the  system's  lifetime  distribution. 
When  we  don't  know  the  lifetime  distribution  explicitly,  the  systems  life  distribution  and 
the  optimal  age  replacement  policy  needs  to  be  estimated.  Estimation  based  on  a  fixed 
sample  of  independent  and  identically  distributed  (i.i.d.)  system  life  times  has  been  ex¬ 
amined  in  detail  [Refs.  5, 6, 7.8].  In  order  to  collect  such  i.i.d.  data,  experimental  systems 
must  be  left  in  service  until  failure.  When  observation  of  complete  system  lifetimes  is 
not  available  (because  it  either  takes  too  much  time  or  is  too  costly),  the  most  cost  ef¬ 
fective  approach  is  to  start  with  an  initial  estimate  of  the  optimal  replacement  age  and 
then  update  this  estimate  after  each  s\  stem  replacement.  After  each  replacement,  the 
next  s\  stem  is  subject  to  a  replacement  polic}  that  is  close  to  the  best  estimated  polic\ 
so  far  [Ref.  9:  pp.  2-3],  This  procedure  is  described  in  detail  in  Chapter  11.  We  will  use 
simulation  to  stud\  this  sequential  estimation  procedure  when  the  under!} ing  life  dis¬ 
tribution  comes  from  a  parametric  family.  In  particular,  in  Chapter  III  we  study  this 
procedure  with  an  underlying  Weibull  distributions  and  in  Chapter  IV  we  consider 
Gamma  life  distributions.  Conclusions  and  recommendations  are  given  in  Chapter  V. 
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II.  THE  SEQUENTIAL  ESTIMATION  PROCEDURE 


A.  THE  OPTIMAL  REPLACEMENT  AGE 

Consider  a  system  or  component  which  is  replaced  upon  failure.  When  the  cost  of 
a  replacement  that  is  planned  in  advance  is  less  than  the  cost  of  an  unplanned  replace¬ 
ment,  a  simple  age  replacement  policy  can  lead  to  considerable  savings.  In  this  type  of 
policy,  an  age  cp  is  specified;  items  that  are  still  functioning  at  that  age  are  replaced 
(these  are  planned  replacements);  items  which  fail  and  are  thus  replaced  prior  to  q>  are 
the  unplanned  replacements. 

Let  .A',,  X2, ...  be  a  sequence  of  independent  and  identically  distributed  (i.i.d.)  positive 
random  variables  with  common  distribution  function  Fe  belonging  to  a  family  of  dis¬ 
tributions  parameterized  by  6.  The  sequence  Xu  X2, ...  represents  the  sequence  of  system 
lifetimes  that  would  be  observable  if  the  systems  were  replaced  at  failure.  Let  C„  C2.  ( 
C,  >  C2 )  be  the  respecti\e  costs  of  planned  replacement  (before  system  failure)  and  un¬ 
planned  replacement  (at  system  failure).  The  observed  durations  between  replacement, 

min(.V,.  <p)  i  =  1,  2 . form  a  renewal  processes  [Ref.  2:  p.  87].  Therefore,  the  long 

run  expected  cost  per  unit  time  with  age  replacement  at  <p  is 

C,  x  Fn{w)  +  C2  x  F0((p) 

CUp)  =  - - - 

El  min  (A},  <p)  ] 

^  C,  x  /•„(<?)  +  C2  x  FM  (2-1) 

_  * 

W  dx 

o 

where  Fe{x)  —  1  —  Fe(x)  is  the  survival  function.  In  equation  (2.1)  the  numerator  is  the 
expected  cost  of  one  replacement  under  the  age  replacement  policy,  and  the  denominator 
is  the  expected  time  between  replacement. 
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Under  reasonable  conditions,  there  exists  a  unique  and  finite  time  q>*  <  oo  where 
C{(p )  attains  a  global  minimum  [Ref.  10:.  pp.  161-168].  For  example,  a  sufficient  con¬ 
dition  for  the  existence  of  <p*  is  that  Fe  has  a  failure  rate  ?.(t)  that  strictly  increases  to 
infinity. 

B.  THE  SEQUENTIAL  ESTIMATION  PROCEDURE 

Finding  the  optimal  planned  replacement  age  <p*  requires  knowledge  of  the  under¬ 
lying  distribution  function  Fg.  If  the  underlying  lifetime  distribution  F$  is  completely 
specified,  then  finding  <p*  can  be  done  either  explicitly  or  numerically  using  equation 
(2.1).  However,  if  the  parameter  6  is  unknown,  we  need  to  estimate  <p*.  If  full  lifetimes 
A',,  X2,  ....  X„  are  available,  then  C(<p)  and  thus  <p*  can  be  estimated  by  replacing  0  in 
equation  (2.1)  with  the  usual  parametric  estimators  of  6.  This  approach  has  the  obvious 
disadvantage  that  it  requires  several  s\  stems  to  operate  until  failure  to  be  observed  be¬ 
fore  a  replacement  policy  is  implemented.  In  addition,  subsequent  data  is  not  used  to 
update  the  policy.  A  more  practical  and  potentially  more  cost  efiectit  e  approach  is  to 
update  estimate  q>*  after  each  replacement  and  implement  the  updated  age  replacement 
policy  (through  estimates  of  <p*)  after  each  replacement. 

Let  <p„  q>2 ,  ...  be  a  sequence  of  estimators  of  q o  where  (f>„  depends  on  the  first  n  re¬ 
placement  ages.  A  procedure  to  compute  the  estimators  q>n  is  developed  as  follows: 

1.  At  the  «th  replacement  observe  the  system  lifetime  Xn  or  whichever  comes  first. 
Let  Zr  =  min(A"„,  </>„_,)  and  Sn  =  I(X„  <  <p„_,)  where  1(A)  is  an  indicator  function  of 
the  set  A.  In  other  words,  if  the  unit  is  replaced  before  failure,  then  the  replace¬ 
ment  time  Z„  =  <pn_,  and  6r  =  1 ,  otherwise  Z„  =  A'  and  <5„  =  0. 

2.  The  data  available  to  estimate  q>n  are  the  pairs  (Z„  S,)  i=  1,  2,  ...,  n  (Throughout, 
without  loss  of  generality  we  take  special  case  as  Z,  =  A',,^,  =  I).  The  maximum 
likelihood  estimator  (MLEj  6„  of  0  is  computed  from  this  right  censored  data. 

3.  <?„  is  then  found  to  minimize  equation  (2. 1 )  with  taking  the  place  of  the  unknown 
parameter  0. 
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The  procedure  is  then  repeated.  The  technical  problem  of  using  this  data  to  estimate  cp 
at  each  stage  is  that  the  pairs  (Z„  8)  i=  1,  2, ....  n  are  clearly  not  i.i.d.  Thus,  the  usual 
methods  for  studying  the  properties  of  the  sequence  of  estimators  {<?„}  from  the  right 
censored  data  are  not  immediately  applicable.  Thus,  we  will  use  simulation  to  study 
both  the  large  and  small  sample  properties  of  this  sequential  estimation  procedure. 

The  replacement  cost  for  the  /th  system  is  C,  if  X,  <  <p,_„  otherwise  the  replacement 
cost  is  C2.  With  this  sequential  estimation  procedure  the  actual  total  replacement  cost 
for  the  first  n  systems  that  are  observed  is  given  in  equation  (2.2) 

n 

Cn  -  Yj  {C>  X  5‘  +  C2  X  (  1  *  )  >’  (2'2) 

(=1 

and  equation  (2.3)  is  the  total  operating  time  for  the  n  systems, 


n  n 

t.-'Z  ™in ■)- E  z'-  (13) 

(=i  ;-i 

When  studying  the  properties  of  this  sequential  estimation  procedure  it  is  important  to 

£ 

sec  if  and  how  fast  the  actual  cost  per  unit  time  -f-  converges  to  optimal  cost  C{q>*) 

‘n 

(Ref.  11]. 

In  this  thesis  we  simulate  the  sequential  estimation  procedure  for  two  parametric 
families  of  distributions  considered:  In  Chapter  III,  F  belongs  to  a  Weibull  family  of 
distributions  with  the  shape  parameter  a  >  1  and  unknown  scale  parameter  (/),  and  in 
Chapter  IV,  F  belongs  to  a  Ga^una  famih  of  distributions  with  shape  parameter  p  >  1 
and  unknown  scale  parameter  {6).  Both  the  Weibull  and  Gamma  distributions  were 
chosen  for  the  simulation  because  they  arc  commonly  used  to  model  system  lifetimes, 
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they  have  increasing  failure  rates  when  their  shape  parameters  >  1,  and  because  esti¬ 
mation  of  the  unknown  parameters  and  minimization  of  the  estimated  cost  function  are 
numerically  tractable. 
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III.  SIMULATION  SETTING  FOR  THE  WEIBULL  DISTRIBUTION 


A.  UNDERLYING  LIFE  DISTRIBUTION 

This  chapter  is  concerned  with  the  estimation  of  the  optimal  replacement  time  when 
it  is  known  that  the  underlying  lifetime  distribution  is  a  member  of  the  two  parameter 
Weibull  family  with  shape  parameter  a  and  scale  parameter  ).  >  0,  where  the  density  is 
given  by 

J[t)  =  a  X  e~^  for  t  >  0.  (3.1) 

The  Weibull  distribution  has  failure  rate 

/(/)  =  a X  (Xtf~x ,  t  >  0 .  (3.2) 

When  a  >  1.0,  the  failure  rate  in  equation  (3.2)  is  strictly  increasing  to  infinity.  Thus, 
for  Weibull  distributions,  with  a  >  1.0,  a  unique  and  finite  optimal  replacement  age  q>* 
exists.  To  compare  the  results  with  [Ref.  9:  p.  12],  the  same  ten  different  Weibull  dis¬ 
tributions  used  in  the  simulation  have  a  values  1.1,  1.2,  ...,  1.9,  2.0.  This  selection  of  a 
values  gives  us  a  range  of  distributions  which  become  more  like  the  exponential  distrib¬ 
ution  as  a  decreases  from  2.0  to  LI.  The  Weibull  distribution  with  shape  parameter 
a  =  2.0  is  called  the  Ra\leigh  distribution  and  has  a  linearly  increasing  failure  rate.  To 
make  fair  comparisons  between  Weibull  distributions,  the  scale  parameter  X  was  chosen 
so  that  the  expected  system  lifetime  E(A')  =  2.0.  In  our  figures  we  have  selected  only 
certain  values  (a  =  1.2,  1.4, 1.6, 1.8,  2.0)  for  the  scale  parameter  a.  This  was  done  to 
make  the  figures  more  readable.  See  Tigurcs  1  and  2,  for  plots  of  the  Weibull  densities 
and  corresponding  failure  rates. 
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1TEIBULL  DENSITY  FUNCTION 
f(t)  =  aXCXt^-^e-W 


Figure  I.  The  Weibull  Density  Function  f(t)  with  E( A’)  =  2.0 
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WEIBULL  DENSITY  FAILURE  RATE  FUNCTION 


X(t)  =  «x(xt)<*-o 


Figure  2.  The  Failure  Rate  of  The  Weibull  Distribution  with  £(,Y,)  =  2.0 


B.  OPTIMAL  REPLACEMENT  TIME 

When  the  distribution  of  the  system  lifetime  is  Weibull,  the  reliability  (survival) 
function  defined  by 


t  >  o ,  ;.  >  0  ,  (3.3) 

i<  0. 

The  long  run  expected  cost  per  unit  time  under  a  simple  age  replacement  policy  with 
scheduled  replacement  at  age  <p,  by  using  equation  (2.1)  is  given  by 


m  =  p(x>t) 
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Q,»  = 


Cj  ( 1  —  r<^),)  +  c2r<^>“ 


r 


,-<**>* 


(3.4) 


</x 


See  Figure  3,  for  plot  of  the  long  run  expected  average  cost  function  when  the  underly¬ 
ing  life  distribution  is  Weibull  with  shape  parameter  a  varying  from  1.1  to  2.0.  For  each 
curve  on  Figure  3  the  optimal  replacement  time  <p*  can  be  located  on  the  x-axis  at  the 
min’mum  point. 


Figure  3.  The  Long  Run  Expected  Average  Cost  Curves  with  £(A',)  =  2.0 
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C.  SEQUENTIAL  ESTIMATION  PROCEDURE  OF  THE  WEIBULL 
DISTRIBUTION 

We  will  use  the  sequential  estimation  procedure  which  we  discussed  in  general  in 
Chapter  II,  with  the  known  shape  parameter  a  >  1.0  and  the  unknown  scale  parameter 


Let  /„  be  the  MLE  of  ).  after  n  replacements.  The  corresponding  estimator  <p„  of 
cp*  which  minimizes  0;  J^cp)  needs  to  be  found  numerically.  The  following  result,  sim¬ 
plifies  this  procedure  considerably. 

Lemma  7.  Let  cp*  minimize  Q  ,(</>)  then  q>*  =  /.x*  where  x*  minimizes  Ci,,(jf). 

PROOF.  By  using  equation  (2.1),  the  cost  function  at  age  t  with  scale  parameter  /  be¬ 
comes, 


Q,*(  0  - 


ri  oo 

J{x)  dx  +  C2  j[x) , 

_ 

rtroc 

A f)  dv  du 

vrt  J.j 


If  we  substitute  the  Weibull  density  into  the  equation  (3.5),  the  numerator  can  then  be 
shown  to  be 


C,  /’W3  1  e~  ^  dx  +  C2  /.^ax*  ]c  ^  dx, 


with  the  change  of  variable  v  =  /.x.  we  obtain 
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(3.6) 


If  we  factor  the  constant  term  (  —  )*”!  from  (3.6),  we  get 


— pr  |ci Jo  ay*~'e~y  dy  +  «/  y  4hJ-  P-7) 

In  expression  (3.7)  the  constant  terms  cancel,  and  the  integrands  are  just  the  density  of 
the  Weibull  distribution  with  /  =  1.  Thus  (3.7)  is 

C,F,=1(;.t)  +  c2f,=1(//).  (3.8) 

With  the  change  of  variables  y  =  /v  the  integrand  in  the  denominator  of  equation  (3.5) 
becomes, 


ay*  ]e~y  dy  =  FM(Xt). 
h.u 

From  (3.8)  and  (3.9),  the  cost  function  may  be  written  as 


c,Jj) 


+  Q^;.=i  (/'-0 

’ i 

FJ=](/u)du 

^0 


(3.9) 


Set  y  =  /u,  then 
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(3.10) 


C,M) 


=  iClm 


A— 1. 


Therefore,  if  x*  minimizes,  CMt,(t)  and  <p*  minimizes  then  <p*  =  Ax*.  ■ 

To  estimate  A.  (unknown  scale  parameter)  from  the  right  censored  data;  recall  that, 

A 

the  available  data  to  estimate  A.„  are  the  pairs  of  (Z„  5.)  i=  1,  2,  ...»  n  where  Z,  =  A'„  and 
<5,  =  0.  By  using  the  maximum  likelihood  estimation  (MLE)  procedure  we  can  estimate 

A 

/,  from  the  first  observation  X\.  The  MLE  of  A.  can  be  found  by  differentiating  the 
likelihood  (or  the  loglikclihood)  with  respect  to  /,  setting  the  results  equal  to  zero,  and 
solving  for  /. 

L(A,  a|Z)  =  o.A.\xxf~xe~^\ 


/(/,  a| Z)  =  In  a  +  a  In  /  +  (a  —  1)  In  xl  —  A.ax*, 


cl{)„  ot\Z) 

CA 


1 

JC 


a 


ei(A  a|Z)  A  A  I 

- - -  =  o  => 

C/.  1  A1 

By  Lemma  1,  we  can  estimate  q>*  such  that  q>*  =  where  <p*  minimizes  G  (<p) 

*!’  * 

and  x*  minimizes  C,  ,(jc). 

After  the  second  observation,  we  set  Z2  —  min(A'2,  <£>,*)  and 

f  1  X2<v* 

L  0  otherwise. 
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After  n  observations,  data  will  be  (Z„  <5,),  (Z2,  S2) .  (Z„,  Sn)  where  Z,  =  min(A'  , 

$*i),  Xt  is  the  /til  lifetime,  and  </>*,  is  the  best  estimation  of  our  optimal  policy  so  far. 


In  general,  if  we  repeat  the  sequential  estimation  procedure  n  times,  we  obtain 
L(/,  «|Z)  =  M)  C^)]'"*2 ...  lf[zn)l6"  [F(z„)]1_<5n, 

/(/,  oc|Z)  =  «(1  +  <52  +  ...  +  8n)  In  ).  —  (z“  +  z2  +  ...  +  z%), 


A 

By  Lemma  1,  wc  can  find  the  MLE  <p„  of  <p*  explicitly  from  the  MLE  of  ).  based 
on  the  first  n  replacements  using 

A  *  0  *  . 

<Pn  =  *nxx  •  (3.12) 
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D.  SIMULATION  RESULTS  FOR  THE  WEIBULL  DISTRIBUflON 


1.  Finding  of  The  Scale  Parameter  ). 

As  mentioned  before,  the  ten  different  Weibull  distributions  used  in  the  simu¬ 
lation  have  a  values  1.1,  1.2,  ....  1.9,  2.0  and  the  scale  parameter  2  was  chosen  so  that 
the  expected  system  lifetimes  £(A')  =  2.0.  The  mean  of  the  Weibull  distribution  is 


£[A']  = 


n±) 

a X 


(3.13) 


In  equation  (3.13),  the  scale  parameter  X  value  can  be  obtained  easily  since  the  scale 
parameter  a  and  £[A'H  are  known.  The  APL  program  Weibull  in  Appendix  A  calculates 
the  scale  parameter  X  values  for  given  a  and  The  results  are  given  in  Table  1  for 

the  shape  parameters  and  corresponding  scale  parameters. 


Table  1.  SCALE  PARAMETER  /  VALUES  FOR  £(AQ  =  2.0 


Shape  pa¬ 
rameter  a 

a  2.0 

a  -  1.9 

a  =  1.8 

a  =  1.7 

a  —  1.6 

Scale  pa¬ 
rameter  X 

0.44311346 

0.44368166 

0.44464337 

0.44612225 

0.44828714 

Shape  pa¬ 
rameter  a 

a  =  1.5 

a  =  1.4 

o»1.3 

a  =  1.2 

m 

Scale  pa¬ 
rameter  X 

0.45137265 

0.45571117 

0.46178836 

0.47032793 

0.48245624 
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2.  Finding  of  The  Actual  Optimal  Age  Replacement  Time  and  Corresponding  Actual 


Cost 

The  actual  age  replacement  time  <p*  can  be  located  in  Figure  3  on  the  x-axis 
at  the  minimum  point  of  the  cost  function  C^,(<p)  (3.4).  In  order  to  find  the  unique 
minimum  point  <p*,  it  hard  to  solve  equation  3.4  analytically.  Therefore,  the  actual  age 
replacement  time  (p *  found  by  simulating  the  cost  function  (3.4)  for  different  cp  where 
the  simulated  cost  function  is  given  by 


C,  x 


|  £/(*,<:  n>) 


r  £/«£„>)  - 

+  C,x  - - ] 


f 


^min(A),  (p) 


(3.14) 


where  X,  (i=  1,  2.  ...,  n)  are  the  simulated  /./.</.  Weibull  random  variables  with  the  spe¬ 
cific  a  and  /.  cp*  is  the  minimum  of  C,, „(</>).  Since  the  optimal  replacement  time  <p* 
comes  from  the  simulation  result,  it  varies  slightly  with  the  number  of  pseudo  random 
variables  used  and  the  seed  numbers  used  to  generate  them,  thus  wc  choose 
n  =  15  x  10s. 

3.  Finding  of  The  Average  Age  Replacement  Times  and  The  Average  Costs 

In  order  to  find  the  average  of  the  estimated  optimal  age  replacement  times  and 
the  average  actual  costs  based  on  lifetimes  with  a  Weibull  distribution,  we  wrote  the 
Fortran  simulation  program  Avcrweib  in  Appendix  C.  Trom  Lemma  1,  wc  know  that 
x*  minimizes  Cj.,<a(/)  and  <p*  minimizes  Q,(/)  when  ip *  =  ).x* .  Thus,  x *  can  be  found 
dividing  the  actual  optimal  age  replacement  times  by  the  scale  parameter  /,  where  the 
actual  age  replacement  time  is  taken  from  the  simulation  program  Sim  and  the  /  values 


are  taken  from  the  Table  1  for  the  specific  value  of  the  shape  parameter  a.  The  un¬ 
planned  and  planned  replacement  costs  C,  and  C2  (to  calculate  the  cost  function);  the 
scale  parameter  ).  (to  calculate  the  x*  value);  the  actual  age  replacement  time  (from 
program  Sim  to  calculate  the  x*  value  and  to  find  MSE  for  each  estimated  age  replace¬ 
ment  time);  the  actual  cost  value  (from  program  Sim  to  find  MSE  value  for  each  esti¬ 
mated  cost);  and  the  shape  parameter  a  must  be  given  by  the  user  in  the  initialization 
part  of  the  simulation  program  Awerweib.  In  much  of  the  simulation,  we  considered  the 
sample  sizes  A*  =  10,  50,  250  and  1000.  Other  sample  sizes  are  also  considered,  but  in 
much  less  detail.  Each  repetition  of  the  simulation  is  based  on  generating  A'  system 
lifetimes.  These  system  lifetimes  (A')  are  used  one  at  a  time  for  the  sequential  estimation 
procedure.  At  the  first  observation,  Z,  =  A'„  and  <5,  =  1.  The  unknown  scale  parameter 

A 

A  values  are  estimated  by  using  the  equation  (3.11).  After  finding  A,  the  estimated  age 
replacement  time  values  are  calculated  by  using  the  equation  (3.12)  in  the  simulation. 
By  using  .V  generated  system  lifetimes,  we  determine  A'  estimated  age  replacement  times 
and  A’  estimated  costs  in  each  simulation.  We  repeat  this  simulation  1000  times 
(NREP=  1000)  and  then  we  find  the  average  value  for  both  age  replacement  times  and 
replacement  costs. 

Let  C„v,  j  =  1,  2,  ...,  1000  be  the  actual  cost  per  unit  time  for  the  first  A’  re¬ 
placements  of  the  jth  repetition  of  a  simulation  (where  each  CjN  is  computed  using  (2.2) 
and  (2.3)).  The  average  actual  cost  per  unit  time  over  the  1000  repetitions  is 


C,v  = 


(3.15) 


IS 


Let  <p*,j  =  1,2, ...,  1000  be  the  estimated  optimal  replacement  time  for  the  first 
N  replacements  of  the  jth  repetition  of  a  simulation.  The  average  age  replacement  time 
over  the  1000  repetitions  is 


1000 


7=1 


(3.16) 


For  each  simulation,  we  also  calculated 


1000  /■», /  ^f\\2 

Z(C,v-C(<p  )) 

^-iooo - • 


7=1 


(3.17) 


M  STAGE 


1000  ,  *  *v2 


Z 

7=1 


1000 


(3.18) 


where  MSECOS  is  the  a\erage  squared  difference  of  the  actual  replacement  cost  per  unit 
time  from  the  estimated  minimum  long  run  expected  replacement  cost  per  unit  time  and 
M STAGE  is  the  average  squared  difference  of  the  actual  optimal  age  replacement  time 
from  the  estimated  optimal  age  replacement  time.  These  MSE  values  are  calculated  in 
the  simulation  in  order  to  see  if  the  sequential  estimation  procedure  is  converging  the 
actual  cost  and  the  actual  optimal  age  replacement  time.  If  we  get  the  MSE  values  close 
to  the  zero,  then  we  can  say  that  this  procedure  is  working  well. 

Tables  2,  3,  4  and  5  summarize  the  simulation  results  of  the  optimal  age  re¬ 
placement  times  and  the  minimum  long  run  expected  replacement  costs  per  unit  time  for 
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different  values  of  a  with  N=  1000  for  fixed  costs  C(  =  2.0,  C,  =  5.0,  C,  =  8.0,  C,  =  10.0 
while  holding  the  C2  fixed  at  1.0.  Tables  6,  7  and  8  summarize  the  simulation  results 
using  the  three  sample  sizes  of  10,  50  and  250  as  small,  moderate  and  large  sample  sizes 
for  different  values  of  a  with  fixed  costs  C,  =  2.0  and  C,=  1.0.  Included  in  Tables 
2~8,  is  the  probability  that  a  system  will  fail  before  the  time  <p  under  the  optimal  age 
replacement  policy, 


P{Xi<<p*)  =  )  .  (3.19) 

We  have  also  plotted  the  results  of  the  average  age  replacement  times  (Fig.  4) 
and  the  average  costs  (Fig.  5)  for  a  -  1.2,  1.4,  1.6,  1.8,  2.0  when  fixed  costs  C,  =  5.0  and 
C2  =  1.0.  From  these  plots  and  the  results  from  Tables  2~S,  we  observe  that  in  general 
when  a  decreases  from  2.0  (i.e.,  the  underlying  life  distribution  is  becomes  more  expo¬ 
nential)  the  optimal  age  replacement  time  <p*,  the  long  run  expected  optimal  replace¬ 
ment  cost  and  MSE  values  increase.  For  values  of  a  close  to  1.0.  very  little  is  gained  by 
replacing  the  system  before  failure  at  the  higher  cost  C,.  When  the  ratio  of  the  un¬ 
planned  replacement  cost  C,  to  the  planned  replacement  cost  C2  increases  the  optimal 
replacement  time  for  a  specific  a  decreases.  The  larger  values  of  (p*,  insure  that  a  small 
percentage  of  replacement  will  be  made  before  failure,  which  is  what  we  desire  if  the 
system's  life  distribution  is  close  to  exponential. 

By  looking  at  the  tables  for  different  sample  sizes,  we  also  observe  that  the  av¬ 
erage  cost  per  unit  time  up  to  the  A’th  replacement  decreases  with  A',  the  number  of  re¬ 
placement.  This  result  is  promising  because  the  ultimate  goal  of  the  sequential 
estimation  procedure  is  to  decrease  costs  while  sampling.  Even  though  as  A'  -*  oo,  the 
long  run  average  cost  per  unit  time  will  approach  the  optimal  replacement  cost  C(<p*) 
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with  probability  1.0  [Ref.  11],  there  is  no  guarantee  that  the  average  replacement  cost 
will  decrease  for  the  first  observations. 
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Table  2. 


ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
WEIBULL  MODEL  WITH  C,  =  2.0,  Q  =  1.0,  E(X)  =  2.0  AND  N=  1000 
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Table  3.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
WEIBULL  MODEL  WITH  C,  =  5.0,  C2  =  1.0,  E( X,)  =  2.0  AND  N=  1000 
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Table  4.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
WEIBULL  MODEL  WITH  C,  =  8.0,  Ct  =  1.0,  E(X.)  =  2.0  AND  N=  1000 
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Table  6.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
WEIBULL  MODEL  WITH  C,  =  2.0,  C2 «  1.0,  AND  £(<¥,)  =  2.0  (N  =  10) 


26 


27 


Table  8.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
WEIBULL  MODEL  WITH  C,  =  2.0.  C,  =  1.0,  AND  £( A',)  =  2.0  (N  =  250) 


2S 


g  run 


Figure  4.  The  Average  Age  Replacement  Times  For  The  Weibull  Model 
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IV.  SIMULATION  SETTING  FOR  THE  GAMMA  DISTRIBUTION 


A.  UNDERLYING  LIFE  DISTRIBUTION 

This  chapter  is  concerned  with  the  estimation  of  the  optimal  replacement  time  when 
it  is  known  that  the  underlying  lifetime  distribution  is  a  member  of  the  two  parameter 
Gamma  family  with  shape  parameter  p  >  1  and  scale  parameter  6  >  0,  where  the  density 
is  given  by 

f[t)  =  — - - — ~~  for  t>  0.  (4.1) 

JKl  T{p)Bp  J 

Up  =  1.0,  the  Gamma  density  with  scale  parameter  6  reduces  to 

f;)=je-7  t>0.  (4.2) 

In  equation  (4.2),  we  have  an  exponential  density  with  parameter  6.  The  reciprocal  of 
the  failure  rate  of  the  Gamma  distribution  is  given  in  equation  (4.3) 

To  =  \fl-r' '~7r  dx,  ,6  0.  (4.3) 

By  solving  equation  (4.3)  analytically  for  different  values  of  p,  we  get  useful  functional 
forms  of  r(t)  in  Table  9  on  page  31.  The  failure  rate  for  the  Gamma  random  variable 
with  p>  1.0  is  a  strictly  increasing  continuous  function  and  is  bounded  above  by  (-^-). 
A  unique  and  finite  optimum  replacement  policy  cp*  exists  and  will  be  finite  if  and  only 
if  [p  —  1)  is  strictly  greater  than  -p:  —  ,  where  C,  and  C2,  are  unplanned  and  planned 
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replacement  costs  (C,  >  C2).  Similar  to  the  Weibull  case,  this  unique  minimum  of  C(<p) 
occurs  at  the  minimum  point  wheic  the  first  derivative  is  zero. 


Table  9.  FUNCTIONAL  FORM  OF  r{t) 


Shape  param¬ 
eter  p 

Scale  parame¬ 
ter  6 

Functional  form 

2.0 

1.0 

t 

t+  1 

3.0 

0.666667 

216 

18/2+  24/+16 

4.0 

0.5 

8/3 

4r3  +  66  +  6t  +3 

5.0 

0.4 

3 1 256 

1310/*  +  2000F  +  2400F  +  1920/  +76S 

In  this  chapter  cost  C,,  C  are  chosen  so  that  <p*is  finite  and  unique.  We  used  shape 
parameter  p  =  2.0.  3.0,  4.0,  5.0  (/>>  1)  and  scale  parameter  6  =  22,  2.3,  2/4.  2,5  in 
our  simulation,  respectively.  This  selection  of  p  values  gives  us  a  range  of  distributions 
which  become  more  like  the  exponential  as  p  decreases  from  5.0  to  2.0.  To  make  fair 
comparisons  between  Gamma  distributions,  the  scale  parameter  6  was  chosen  so  that 
the  expected  system  lifetime  L(A')  =  2.0.  See  Figures  6  and  7,  for  plots  of  the  Gamma 
densities  and  corresponding  failure  rates. 
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GAMMA  DENSITY  FUNCTION 
f(t)  =  /  r(p)9» 


Figure  6.  The  Gamma  Density  Function  f(t)  with  E{ Xf)  =  2.0 
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GAMMA  DENSITY  FAILURE  RATE  FUNCTION 


Figure  7.  The  Failure  Rate  of  The  Gamma  Distribution  with  £(A',)  =  2.0 

B.  OPTIMAL  REPLACEMENT  TIME 

When  the  distribution  of  a  system  lifetime  is  Gamma,  the  expected  long  run  a\era'ge 
cost  per  unit  time  under  a  simple  age  replacement  policy  with  scheduled  replacement  at 
age  t.  b\  using  equation  (2.1)  is  given  by 


G 


<*<  _  . 
xp  1  e  o 


n p)0P 


dx  +  C2 


r°°  n_, 

.r  e  o 
V{p)0p 


dx 


(4.4) 


l\x)  dx 
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In  equation  (4.4),  the  corresponding  survival  function  F(t )  is, 


F(/)  = 


„  1  * 
Xs?  ]e  e 


Tip)^ 


dx, 


i  >  0. 


(4.5) 


In  this  thesis,  we  obtained  the  numerical  values  for  F{i),  by  solving  equation  (4.5)  ana¬ 
lytically  for  some  specific  values  of  p.  Table  10  gives  the  useful  functional  form  of 
F{t).  See  Figure  8  for  plot  of  the  Qt/>(/)  when  the  underlying  life  distribution  is  Gamma 
with  shape  parameter  p  varying  from  2.0  to  5.0.  For  each  curve  on  Figure  8  the  optimal 
replacement  time  <p *  can  be  located  on  the  x-axis  at  the  minimum  point. 


Table  10.  FUNCTIONAL  FORM  OF  F(T) 


Shape  param¬ 
eter  p 

Scale  parame¬ 
ter  0 

Functional  form 

2.0 

1.0 

t  +  1 

<?f 

3.0 

0.666667 

2.25F  4-  3/  4-  2 

2<?l~' 

4.0 

0.5 

4/5  +  6/J  4-  6/  +  3 

3e20: 

5.0 

0.4 

1310/4  +  2000/3  +  2400/J  +  1920/  4-  768 

768<?:~r 
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LONG  RUN  EXPECTED  AVERAGE  COSTS  FOR  GAMMA  DISTRIBUTION 


(Cj  =  5.0  ,  C,  =  1.0) 


Figure  8.  The  Long  Run  Expected  Average  Cost  Curves  with  £(A')  =  2.0 


C.  SEQUENTIAL  ESTIMATION  PROCEDURE  OF  THE  GAMMA 
DISTRIBUTION. 

For  the  Gamma  distribution,  we  will  use  our  parametric  sequential  estimation  pro¬ 
cedure  again  which  we  mentioned  in  Chapter  II,  with  the  known  shape  parameter 
p  >  1.0  and  the  unknown  scale  parameter  0. 

A 

As  in  the  Weibull  case,  minimizing  G  (/)  (where  0n  is  the  MLE  of  6  after  n  re- 

^  *n'P 

placements)  can  be  simplified  with  the  following  result: 

X* 

Lemma  2,  Let  <p*  minimize  Q,(<p)  then  <p*  =  —r~  where  x*  minimizes  Cl-P(.v). 
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PROOF.  From  the  cost  function  CBtP(<f)  for  the  Gamma  distribution  in  equation  (4.4), 
after  making  the  substitution  y  =  -r,  the  numerator  may  be  written  as 


e  yP-'QPe-y 
1  OTip)^ 


d  dy  + 


/*  oo  n  _ 

2  er(p )dp 

%)  i 


0  dy. 


Canceling  the  scale  parameters  (6)  we  obtain 


(4.6) 


C, 


/»— 

T'e-* 


f*oc 


^0 


!"(/>) 


dy  +  C2 


V' 


n») 


dv. 


(4.7) 


The  integrands  of  equation  (4.7)  are  the  density  of  the  Gamma  distribution  with  6=1. 
Finally  the  numerator  of  (4.4)  becomes 


C./Wy)  +  C2Fd=I(f).  (4.8) 

In  equation  (4.4)  the  denominator  may  be  written  as 


dv  du. 


(4.9) 


If  we  first  solve  the  inner  integral  b>  changing  variable  y  with  y,  vve  get 
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[  J{v)dv  = 


OO  v 

/“Vi 


r(p)^ 


p-i 


T(p)^ 


After  cancellation,  the  denominator  is 


r—  >  —v 

r  e  ' 


(4.10) 


Again  the  integrand  of  equation  (4.10)  is  the  density  of  the  Gamma  distribution  with 
0=1,  From  (4.8)  and  (4.10),  the  cost  function  may  be  written  as 


CgJj)  = 


QW  J)  +  c2FUj ) 


Wf )'« 


(4.11) 


If  we  let  the  y  =  — ,  equation  (4.11)  becomes, 
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Therefore,  if  jc*  minimizes,  Ctf=li,(/)  and  (p *  minimizes  Q,(/)  then  (p*  =  —^~.  ■ 

There  is  no  closed  form  solution  for  the  MLE  of  6  based  on  right  censored  data. 
Thus,  we  use  the  EM  algorithm  to  find  the  MLE  numerically. 

1.  EM  Algorithm 

"  The  Expeciaiion-Maximizaiion  (EM)  algorithm  is  an  iterative  algorithm  used 
to  compute  the  MLE  in  incomplete-data  problems.  The  algorithm  is  applicable  to  more 
general  missing-data  patterns,  but  usually  involves  more  computations  than  methods 
designed  for  special  patterns.  The  iterations  of  the  EM  algorithm  are 

1.  Replace  missing  values  by  estimated  values. 

2.  Estimate  parameters. 

3.  Recstimatc  the  missing  values  assuming  the  new  parameter  estimates  are  correct. 

4.  Recstimatc  parameters. 

and  so  forth,  iterating  until  convergence  "  [Ref.  12:  pp.  127-141).  Since  each  iteration 
of  the  algorithm  consist  of  an  expectation  step  (E)  followed  by  a  maximization  step  (M) 
it  calls  the  EM  algorithm. 

The  M  step  is  performed  by  maximizing  the  likelihood  as  if  there  were  no  miss¬ 
ing  data.  Thus,  the  M  step  of  EM  uses  the  identical  computational  methods  as  MLE 
from  1(0,  p\X)  [Ref.  12:  pp.  127-141).  With  the  shape  parameter  p  assumed  known,  the 
maximum  likelihood  estimator  of  0  based  on  n  unccnsored  observations  ,V,.  ....  Xn  is 


(4.13) 


The  E  step  finds  the  conditional  expectation  of  the  missing  (censored)  data  given 
the  observed  data  and  current  estimated  parameters,  and  then  substitutes  these  expec¬ 
tations  for  the  missing  (censored)  data  [Ref.  12:  pp.  127-141],  Thus  at  the  i'th  iteration 

A 

of  the  EM  algorithm  the  MLE  6„  based  on  n  replacements  is  approximated  by 


A 

Ki  = 


n 

YjldiZi  +  (\-6i)E(X\X>Zt)'] 


(4.14) 


where  E{X\X  >  Zj  is  the  conditional  expectation  of  the  random  variable  X  ~  Gamma( 
0  =  The  functional  forms  of  this  conditional  expectation  for  the  specific  shape 

parameter  p  are  given  in  Table  1 1  on  page  40. 

As  in  the  Weibull  case,  we  have  no  information  about  the  age  replacement  time 
(p  at  the  first  observation.  Thus,  the  Z,  value  equals  the  first  observation  A',.  We  find 
0 1  (by  using  equation  (4.13)).  we  get 


•>-  -f- 


From  the  Lemma  2,  we  can  estimate  rp*  such  that  <p*  =  — ,  where  <p*  minimizes 


Cj  ^(<p).  and  .v*  minimizes  Cli(,t.v). 

At  the  second  observation,  we  have  two  cases.  If  Aj  is  less  than  <p*,  then  02  is 
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Table  11.  FUNCTIONAL  FORM  OF  E(X|X><?)  WHERE  A’  ~  GAMfKp) 


Shape  pa¬ 
rameter  p 


Scale  pa¬ 
rameter  6 


2.0 


1.0 


3.0 


0.666667 


4.0 


0.5 


5.0 


0.4 


<r 


Functional  Form 


cp2  +  26<p  +  2  62 
(p  +  0 


(£_  +  3  fl(y2  +  2  Ocp  +  2d1) 
(p2  +  29(p  +  262 


4>4  +  4  0U£  +  3fl(<?2  +  2  6<p  +  2fl2)) 
4 p 3  +  3  0(<p2  +  20  q>  +  202) 


+  50 (<p4  +  ABjcf  +  3fl(y2  +  2  0cp  +  2ft2))) 
4>4  +  4ft(4>3  +  3ft(c>2  +  26ip  +  202)) 


02 


X)  +  A  2 
2/> 


Otherwise,  the  observation  is  censored,  and  we  use  the  equation  (4.14)  for  the  E  (ex¬ 
pectation)  step  of  the  EM  algorithm  until  02  convergence.  These  iterations  arc 


$2,  i  ~ 


0 2,2  ~ 


x,  +  e,U2ix2>^] 

2 /> 


Z,  +  E'2U’2IA2>4>n 
2/> 


and  so  on.  Finally,  these  iterations  converge  02„.  When  the  difference  of  the  absolute 

A  A 

value  of  0 2„,  and  02jn_,  is  small,  the  stopping  criteria  is  satisfied.  At  that  point,  we  can 

A  A  A  A  V* 

replace  the  02„  with  02.  Again,  we  can  estimate  q>*  such  that  q> *  =  -7—. 

The  procedure  is  then  repeated.  After  determining  6  values  for  each  replace¬ 
ment,  we  can  apply  this  estimated  6  values  to  the  Lemma  2.  The  estimated  optimal  age 
replacement  time  can  be  expressed  by, 


■£.  (4-15) 

On 

For  large  «,  the  estimated  optimal  age  replacement  time  <p*  converges  to  an  optimal  age 
replacement  time  40*. 

D.  SIMULATION  RESULTS  FOR  THE  GAMMA  DISTRIBUTION 

1.  Finding  of  The  Scale  Parameter  0 

The  four  different  Gamma  distributions  used  in  the  simulation  have  p  values  2.0, 
3.0,  4.0.  5.0,  and  the  scale  parameter  9  was  chosen  so  that  the  expected  system  lifetimes 
E(A’)=  2.0.  The  mean  of  the  Gamma  distribution  is  E(X)  =  p6.  The  0  value  can  be  ob¬ 
tain  since  both  E(X)  and  p  are  known.  For  example,  if  p-  3.0,  then  6  =  0.666667. 

2.  Finding  of  The  Actual  Optimal  Age  Replacement  Time  and  Corresponding  Actual 
Cost. 

Similar  to  the  Weibull  case,  we  found  the  minimum  value  of  the  cost  function 
(4.4)  and  corresponding  actual  age  replacement  time  q>*  for  different  costs  C,,  C2  and 
parameters  ( 0,p )  by  simulation. 

3.  Finding  of  The  Average  Age  Replacement  Times  and  The  Average  Costs 

As  in  the  Weibull  case,  in  order  to  find  the  average  of  the  estimated  optimal  age 
replacement  times  and  the  average  actual  costs  based  on  lifetimes  with  a  Gamma  dis- 
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tribution,  we  wrote  the  Fortran  simulation  program  Avergam  in  Appendix  D.  From 
Lemma  2,  x*  =  (p*d  where  the  actual  age  replacement  time  cp*  is  taken  from  the  simu¬ 
lation  program  Sim  and  the  6  values  are  calculated  from  the  equation  E (X)=pd  for  the 
specific  value  of  the  shape  parameter  p.  The  unplanned  and  planned  replacement  costs 
C,  and  C2  (to  calculate  the  cost  function);  the  scale  parameter  8  (to  calculate  the  x * 
value);  the  actual  age  replacement  time  (from  program  Sim  to  calculate  the  x*  value  and 
to  find  MSE  for  each  estimated  age  replacement  time);  the  actual  cost  value  (from  pro¬ 
gram  Sim  to  find  MSE  value  for  each  estimated  cost);  and  the  shape  parameter  p  must 
be  given  by  the  user  in  the  initialization  part  of  the  simulation  program.  In  much  of  the 
simulation,  we  considered  the  sample  sizes  A*  =  10,  50,  250  and  1000.  The  other  sample 
sizes  are  also  considered,  but  in  much  less  detail.  Each  simulation  is  based  on  generating 
1000  sequences  of  A  sxstem  lifetimes.  We  used  these  system  lifetimes  (A’)  one  at  a  time 
for  the  sequential  estimation  procedure.  At  the  first  observation  Z,  =  A„  and  d,  =  1. 
As  mentioned  before,  the  scale  parameter  0  values  are  estimated  by  the  EM  algorithm. 
After  finding  6,  the  estimated  age  replacement  time  values  are  calculated  by  using  the 
equation  (4.15)  in  the  simulation.  Similar  to  the  Weibull  case,  by  using  A'  generated 
system  lifetimes,  we  determine  A'  estimated  age  replacement  times  and  A'  estimated  costs 
in  each  simulation.  We  repeat  this  simulation  1000  times  (NREP=  1000)  and  then  we 
find  the  axerage  \aluc  for  both  age  replacement  times  and  replacement  costs.  See  the 
equations  (3.15)  and  (3.16)  of  Chapter  III.  At  each  simulation  we  also  calculated  the 
MSE  values  for  both  the  age  replacement  times  and  the  long  run  expected  average  costs 
by  using  the  equations  (3.17)  and  (3. IS). 

Tables  12,  13,  14  and  15  summarize  the  simulation  results  of  the  optimal  age 
replacement  times  and  the  minimum  long  run  expected  replacement  costs  per  unit  time 
for  different  \aiues  of  p  with  A’=  10O0  for  fixed  costs  C,  =  2.0,  C,  =  5.0,  C,  =  8.0, 
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C,  =  10.0  while  holding  the  C2  fixed  at  1.0.  Tables  16,  17  and  18  summarize  the  simu¬ 
lation  results  using  the  three  sample  sizes  of  10,  50  and  250  with  different  values  of  the 
shape  parameter  p  for  fixed  costs  C,  =  2.0  and  C2  =  1.0.  Included  in  Tables  12~  1 8,  is 
the  probability  that  a  system  will  fail  before  the  time  cp  under  the  optimal  age  replace¬ 
ment  policy  (3.19). 

We  have  also  plotted  the  results  of  the  average  age  replacement  times  and  the 
average  costs  for  different  shape  parameter  p.  See  Figure  9,  for  plot  of  the  average  age 
replacement  times  for  =  3.0,  4.0,  5.0  (p-  2.0  was  not  selected  because  its  average  age 
replacement  times  were  high  acording  to  the  others)  and  Figure  10,  for  the  average  costs 
for  p=  2.0,  3.0,  4.0,  5.0  when  the  fixed  costs  C,  =  2.0  and  C2=  1.0. 

From  these  plots  and  the  results  from  Tables  12 —  1 8,  we  observe  that  when  p 
decreases  from  5.0  to  2.0  (i.e.,  the  underlying  life  distribution  is  becomes  more  expo¬ 
nential)  the  long  run  expected  optimal  replacement  costs  increase.  When  the  ratio  of  the 
unplanned  replacement  cost  C,  to  the  planned  replacement  cost  C2  increases  the  optimal 
replacement  time  for  a  specific  p  decreases.  The  larger  values  of  <p*,  insure  that  a  small 
percentage  of  replacement  will  be  made  before  failure,  which  is  what  we  desire  if  the 
system's  life  distribution  is  close  to  exponential. 

As  in  the  Wcibull  case,  by  looking  the  tables  for  different  sample  sizes,  we  ob¬ 
serve  that  the  average  cost  per  unit  time  up  to  the  A’th  replacement  decreases  with  .V, 
the  number  of  replacement.  Thi  result  is  promising  because  the  ultimate  goal  of  the 
sequential  estimation  procedure  is  to  decrease  costs  while  sampling.  Even  though  as 
.V  -+  oo,  the  long  run  average  cost  per  unit  time  will  approach  the  optimal  replacement 
cost  C(<p*)  with  probability  1.0  [Ref.  11],  there  is  no  guarantee  that  the  average  re¬ 
placement  cost  will  decrease  for  the  first  observations.  For  large  sample  sizes,  our  esti¬ 
mated  average  cp *  and  the  average  C(e*)  values  are  too  close  to  their  actual  values. 
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Table  12.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
GAMMA  MODEL  WITH  C,  =  2.0,  Q  =  1.0,  AND  £(A’,)  =  2.0  (N  =  1000) 


Shape  parameter  p 

p=  2.0 

II 

p  =  4.0 

p=  5.0 

Scale  parameter  6 

1.0 

0.666667 

0.5 

0.4 

x*  =  0x<p* 

35.35100 

2.11607 

1.10275 

0.75876 

Optimal  replacement  time 
<P* 

35.35100 

3.17410 

2.20550 

1.89690 

Average  <p* 

35.37980 

3.17659 

2.20480 

1.89694 

MSE  of  (p* 

0.62197 

0.00371 

0.0016S 

0.00106 

Long  run  expected  optimal 
replacement  cost  C(<p*) 

1.00000 

0.99456 

0.97167 

0  93374 

Average  C[(p*) 

1.000S1 

0.99552 

0.97201 

0.94614 

MSE  of  C(<p*) 

O.OU049 

0.00036 

0.00032 

0.00041 

Pi  A’  <  tp*) 

1.00000 

0.61467 

0.18173 

0.03994 
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Table  13.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
GAMMA  MODEL  WITH  C,  =  5.0.  C  =  1.0.  AND  £U',)  =  2.0  (N=  1000) 


Shape  parameter  p 

■SOI 

XI 

II 

CO 

O 

p  =  4.0 

p=  5.0 

Scale  parameter  0 

1.0 

0.666667 

0.5 

0.4 

x*  =  0xcp* 

1.30500 

0.67220 

0.48920 

0.39952 

Optimal  replacement  time 

<p* 

1.30500 

1.00830 

0.97840 

0.99880 

Average  <p* 

1.30421 

1.00727 

0.97494 

0.99402 

MSE  of  </>* 

0.00177 

0.00107 

0.00106 

0.00202 

Long  run  expected  optimal 
replacement  cost  C(<p*) 

2.26480 

1.87695 

1.63237 

1.44237 

Average  C(</>*) 

2.26S30 

1.8S5S1 

1.64236 

1.48797 

MSE  of  C( <,■)*) 

0.00528 

0.00367 

0.O0307 

0.0 1029 

p>a;<^) 

0.37495 

0.19428 

0.13517 

0.00291 

Table  14.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
GAMMA  MODEL  WITH  C,  =  8.0,  C2  =  l.Q,  AND  £(*,)  =  2.0  (N=  1000) 


Shape  parameter  p 

p=  2.0 

p=  3.0 

II 

o 

p=  5.0 

Scale  parameter  8 

1.0 

0.666667 

0.5 

0.4 

x*  =  0  x  q>* 

0.81890 

0.49540 

0.38620 

0.32956 

Optimal  replacement  time 
<P* 

0.81890 

0.74310 

0.77240 

0.82390 

A'erage  <p* 

0.81677 

0.73945 

0.76560 

0.81510 

MSE  of  (p* 

0.00124 

0.00112 

0.00172 

0.00334 

Long  run  expected  optimal 
replacement  cost  C{(p*) 

3.15166 

2.28405 

1.97650 

1.68813 

Average  C{<p*) 

3.16146 

2.102S1 

2.00406 

1.77106 

MSE  of  C(<p*) 

0.01767 

0.01076 

0.01375 

P(A’  <(p*) 

■i 

0.03960 

0.00SO6 

0.00120 
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Table  15.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 


GAMMA  MODEL  WITH  C,  =  10.0,  C2  =  1.0,  AND  £(A,)  =  2.0 
(N=  1000) 


Shape  parameter  p 

O 

04 

II 

cx 

p=3.0 

T3 

II 

O 

p  =  5.0 

Scale  parameter  0 

1.0 

0.666667 

0.5 

0.4 

x*  =  6  x  <p* 

0.68010 

0.43693 

0.35020 

0.30424 

Optimal  replacement  time 

<p* 

0.68010 

0.65540 

0.70040 

0.76060 

Average  <p* 

0.67754 

0.65103 

0.68916 

0.74911 

MSE  of  </>* 

0.00113 

0.00111 

| 

0.00263 

i 

0.00359 

Long  run  expected  optimal 
replacement  cost  C(<p*) 

3.64327 

2.64540 

2.14725 

1.80642 

Average  C(^*) 

3.66032 

2.67401 

2.19022 

1.90982 

MSE  of  C(y*) 

0.02993 

0.01724 

0.02745 

0.04275 

P(A'  <  q>*) 


892 


0.02894 


0.00576 


0.00082 


Table  16.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 


GAMMA  MODEL  WITH  C,  =  2.0,  C2  =  1.0,  AND  £(A,)  =  2.0  (N=  10) 


Shape  parameter  p 

p=  2.0 

p=  3.0 

II 

o 

II 

Vi 

O 

Scale  parameter  6 

1.0 

0.666667 

0.5 

0.4 

x*  =  0x<p* 

35.35100 

2.11607 

1.10275 

0.75876 

Optimal  replacement  time 
<P* 

35.35100 

3.17410 

2.20550 

1.89690 

Average  <p * 

37.57612 

3.23432 

2.18015 

1.82749 

MSE  of  </>* 

93.99036 

0.51523 

0.24643 

0.16467 

Long  run  expected  optimal 
replacement  cost  Q(p*) 

1.00000 

0.99456  i 

0.97167 

0.93374 

Average 

1.06294 

1.02846 

1.00494 

0.97408 

MSE  of  C((p*) 

0.07521 

0.04433 

0.03426 

0.02609 

P(A :,<(?*) 

1.00000 

0.61467 

0.18173 

0.03994 
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Table  17.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
GAMMA  MODEL  WITH  C,  =  2.0.  C2  =  1.0.  AND  £( Xf)  =  2.0  (N  =  50) 


Shape  parameter  p 

p=  2.0 

HEEHHI 

x> 

II 

o 

p=  5.0 

Scale  parameter  6 

1.0 

0.666667 

0.5 

0.4 

x*  =  dx<p* 

35.35100 

2.11607 

1.10275 

0.75876 

Optimal  replacement  time 

35.35100 

3.17410 

2.20550 

1.89690 

Average  <p* 

35.84104 

3.19201 

2.20190 

1.87172 

MSE  of  <p* 

13.43787 

0.07847 

0.03622 

0.03245 

Long  run  expected  optimal 
replacement  cost  C\<p*) 

1.00000 

0.99456 

0.97167 

0.93374 

Average  C{<p*) 

1.01386 

1.00393 

0.98294 

0.95588 

\ 

MSE  of  C{(p*) 

0.01075 

0.00731 

0.00583 

0.00545 

PC  A',  <v*) 

\ 

1.00000 

0.61467 

0.18173 

0.03994 
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Table  18.  ESTIMATED  OPTIMAL  REPLACEMENT  TIMES  OF  THE 
_ GAMMA  MODEL  WITH  C,  =  2.0,  C\  = 1.0,  AND  £(A~,)  =  2.0  (N  =  250) 


Shape  parameter  p 

p=  2.0 

p=  3.0 

II 

O 

p=  5.0 

Scale  parameter  0 

1.0 

0.666667 

0.5 

0.4 

x*  =  6x(p* 

35.35100 

2.11607 

1.10275 

0.75876 

Optimal  replacement  time 
<P* 

35.35100 

3.17410 

2.20550 

1.89690 

Average  <p* 

35.471 4S 

3.18530 

2.20761 

1.89861 

MSE  of  cp * 

2.53150 

0.01452 

0.00649 

0.00436 

Long  run  expected  optimal 
replacement  cost  C(<p*) 

1.00000 

0.99456 

0.97167 

0.93374 

Average  C(<p*) 

1.00340 

0.99883 

0.97499 

0.94996 

MSE  of  C(<p*) 

0.00202 

0.00141 

0.00124 

0.00127 

P(A’  <<?*) 

1.00000 

0.61467 

0.18173 

0.03994 
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COMPARISON  OF  AGE  REPLACEMENT  TIME 


Figure  9.  The  Average  Age  Replacement  Times  For  The  Gamma  Model 


COMPARISON  OF  AVERAGE  COSTS 
C,=2.0  .  Cg=1.0 


Figure  10.  The  Average  Costs  For  The  Gamma  Model 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


Throughout  this  thesis,  we  have  considered  age  replacement  policies.  In  such  poli¬ 
cies  a  unit  is  always  replaced  at  the  time  of  failure  or  <p  units  of  time  after  its  installation, 
whichever  comes  first.  The  time  q>  is  called  the  scheduled  replacement  time  or  the 
scheduled  censoring  time.  The  optimal  replacement  time  (p*  achieves  the  minimum  long 
run  expected  maintenance  cost  per  unit  time.  The  results  of  the  simulations  show  that 
the  parametric  estimators  work  well  under  the  conditions  for  which  they  were  intended 
and  by  using  the  sequential  estimation  procedure  significant  cost  and  time  savings  can 
be  effected. 

In  most  cases,  the  estimated  optimal  age  replacement  times  and  the  actual  costs  got 
close  to  the  true  optimal  age  replacement  time  and  the  minimum  expected  cost  per  unit 
time  respectively,  even  with  moderate  sample  sizes.  Therefore,  we  conclude  that  our 
sequential  estimation  procedure  of  the  age  replacement  policies  to  minimize  the  long  run 
expected  cost  can  be  applied  to  the  real  problem. 

An  important  part  of  our  analysis,  which  would  be  very  difficult  to  prove  theore¬ 
tically.  showed  that  the  average  actual  cost  per  unit  time  decreased  monotonically  as  the 
sample  si/c  increased.  Such  a  decrease  makes  intuitive  sense,  i.e.  an  age  replacement 
policy  using  an  estimated  <p*  should  do  better  as  more  data  is  collected.  This  result  is 
desirable  for  a  sequential  estimation  procedure. 

Directions  for  Future  Research 

It  is  our  hope  the  ideas  and  the  sequential  estimation  procedure  described  in  this 
thesis  represent  a  plateau  for  the  development  of  the  optimum  age  replacement  policies 
to  minimize  the  long  run  expected  maintenance  costs.  The  minimization  problem  is  in 


principle  difficult  and  future  developments  will  exploit  fundamentally  new  ideas.  Here, 
we  briefly  describe  directions  in  which  future  research  might  be  pursued. 

•  We  have  described  the  scenario  in  which  the  underlying  lifetimes  are  i.i.d.,  along 
with  a  straight  forward  age  replacement  policy.  In  many  situations  this  model  is 
hot  adequate.  For  example,  if  an  item  is  repaired  at  failure  rather  than  replaced, 
the  i.i.d.  assumption  is  equivalent  to  requiring  that  the  repaired  item  function  as 
well  as  a  new  one.  Clearly,  modeling  times  between  failure  of  a  repairable  system 
as  i.i.d.  is  inappropriate.  More  realistic  models  incorporate  the  possibility  that  re¬ 
pairs  are  less  than  perfect.  It  is  also  possible  that  the  quality  of  planned  mainte¬ 
nance  varies  from  time  to  time.  If  an  imperfect  repair  model  is  permitted  how 
should  the  sequential  estimation  procedure  be  changed  to  fit  new  situation. 

•  Sequential  estimation  when  the  underlying  life  distribution  F  comes  from  the 
normal  or  the  modified  extreme  value  distributions  which  have  increasing  failure 
rate. 

•  The  minimization  of  the  long  run  expected  costs  is  not  appropriate  under  all  cir¬ 
cumstances  because  in  this  model  the  planned  and  unplanned  costs  are  fixed. 
Other  cost  functions  need  to  be  considered.  For  example,  costs  can  be  modeled  to 
change  with  time. 

•  How  to  change  our  estimation  procedure  if  we  have  different  systems  such  that 
serial,  parallel  or  bridge  systems  with  same  or  different  lifetimes  rather  than  one 
unit  and  one  lifetime. 
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APPENDIX  A.  APL  CODE  TO  CALCULATE  LAMBDA  OR  MEAN  FOR 

THE  WEIBULL  DISTRIBUTION 


V  WEIBULLiAiEiLiX 

[13  A  THIS  APL  PRO GR AVI  CALCULATES  THE  SCALE  PARAMETER 
[23  A  LAMBDA  AND  THE  E[X3  FOR  THE  WEIBULL  DISTRIBUTION. 
[33  A  THE  SHAPE  PARAMETER  a  AND  THE  OTHER  UNKNOWN 
[43  A  MUST  BE  GIVEN  BY  THE  USER. 

[53  LUCK+0 

[63  START: 1  PLEASE  TYPE  1  IF  YOU  WISH  TO  FIND  LAMBDA ' 

[73  • PLEASE  TYPE  2  IF  YOU  WISH  TO  FIND  £[X3 ‘ 

[83  TYPEl+ft 

[93  +((TYPEl*'l' )*(TYPE1*'2' ))/WARNINGl 

[103  •Kl+ZTP£1=  '  2  '  )/EXPECVA 

[113  LAMBDA: 'PLEASE  ENTER  THE  a  VALUE ' 

[123 

[133  'PLEASE  ENTER  THE  EXPECTED  VALUE  EIX 3 ' 

[143  £*□ 

[153  a  LAMBDA  =  GAMMA (1/a)  /(ax  EIX 3) 

[163  X«-l *A 

[173  C(X-l))*Ux£) 

[183  +END 

[193  EXPECVA:' PLEASE  ENTER  THE  a  VALUE' 

[203 

[213  'PLEASE  ENTER  THE  LAMBDA  VALUE' 

[223  L+ 0 

[233  A  £[X3  =  GAMMA (1/ a)  /(ax  LAMBDA ) 

[243  X+1*A 

[253  (I  (X-l))*(AxL) 

[263  +END 
[273  WARNINGl: 
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[28]  LUCK+LUCK+ 1 

[29]  +(LUCKZ2)/CHEW1 

[30]  ' YOU  CAN  ONLY  ENTER  THE  NUMBERS  1  OR  2  * 

[31]  +START 

[32]  CHEWli x JUST  TYRE  1  OR  2' 

[33]  +START 

[34]  WARNING 2i 

[35]  LUCK+LUCK+1 

[36]  +  UUCKZ2)/CHEW2 

[37]  ’YOP  CM  C/VLY  PWPPP  TAP  LETTERS  Y  OR  N' 

[38] 

[39]  CHEW 2 : 'JUST  TYPE  Y  OR  N' 

[uo] 

[41]  END l 'DO  YOU  WANT  TO  RUN  THE  PROGRAM  AGAIN? (Y/N) * 

[42]  PYPP2+D 

[43]  *((l+YYPE2)='y' )/START 

[44]  -+(1\TYPE2*  1  Y*  )t\(l+TYPE2*N)/WARNING2 

[45]  1  BYE  HYP' 

V 


APPENDIX  B.  APL  CC  ^E  TO  CALCULATE  THE  ACTUAL  AGE 
REPLACEMENT  TIME  AND  CORRESPONDING  THE  ACTUAL 

MINIMUM  COST 


V  SIMiCliC2iIiJiFX-tXiTiXAiYAiC;DiXMINiYMIN 

[I]  fl  THIS  PROGRAM  SIMULATES  THE  COST  FUNCTION  ( EQ .  2.1)  TO 

C 2 ]  R  FIND  THE  MINIMUM  VALUE  ( YMIN )  OF  THE  COST  FUNCTION  AND 

[3]  P  CORRESPONDING  AGE  REPLACEMENT  TIME  ( XMIN )  FOR  THAT 

[4]  p  POINT.  AFTER  FINDING  MINIMUM  VALUES  INSIDE  THE  L00P1 

[5]  fl  IT  REPEATS  THE  PROCEDURE  300  TIMES  INSIDE  THE  L00P2. 

[6]  P  FINALLY,  THE  PROGRAM  GIVES  THE  AVERAGE  VALUES  FOR 

[7]  R  BOTH  MINIMUM  POINT  AS  AXST  AND  ACST . 

[8]  R 

[9]  T+( 15000)4100 

[10:  fl  THIS  GIVES  US  A  VECTOR  OF  T{ 0.01,  0.02,  ....  50) 

[II]  fl  TO  CALCULATE  FIRST  (7(0.01)  AND  THEN  £7(0.02)  UP  TO 
[123  fl  (7(50)  OF  5000  COST  VECTOR. 

[13]  R  INITIALIZATION .. . 

[14]  R  UNPLANNED  AND  PLANNED  REPLACEMENT  COSTS  MUST  BE 

[15]  fl  GIVEN  BY  THE  USER. 

[16]  (71+5 

[17]  (72+1 

[18]  X4+i0 

[19]  I^+i0 

[20]  «7+  0 

[21]  fl  J  IS  THE  INCREMENT  OF  THE  LOOP2  J= 1,  2,  ...»  300 

[22]  fl  MODEL. . . 

[23]  LOOP2: 

[24]  X+5000  WEIRAND  2  2.2567587 

[25]  fl  PREVIOUS  LINE,  GENERATES  5000  SYSTEM  LIFETIMES  FROM 

[26]  fl  WEI(ALPHA=2. 0 ,BETA=2. 2567587)  AS  VECTOR  X.HERE  BETA 
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[27]  r  VALUE  REPRESENTS  1  OVER  LAMBDA-1*0. 44311346 

[28]  A  FOR  GAMMA  DISTRIBUTION  LINE  24  CAN  BE  SNITCH  WITH 

[29]  R  X+5000  GAMRAND  4  0.5  GAM(  P=4  , THETA- 0.5  ). 

[30]  t7+«7+l 

[31]  OlO 

[32]  I<-0 

[33]  R  I  IS  THE  INCREMENT  OF  THE  INNER  LOOP  1=1,  2,  ...,  5000 

[34]  LOOPl; 

[3  5]  1*1+1 

[36]  fl  C  IS  THE  SIMULATED  COST  FUNCTION 

[37]  fl  D  IS  THE  DENUMERATOR  OF  THE  COST  FUNCTION 

[38]  £)•«■(  (+/(XL2’[J]  ))*5000) 

[39]  C*C,  (((C2x(l-FX))  +  (Clx(FX«-((+/X£Z’[I]  >*5000  )  )  )  )*D) 

[40]  A  IN  THE  FIRST  LOOP  C  VECTORS  OBTAIN  FOR  EACH  T 

[41]  ->(J<5000  )/LOOPl 

[42]  YMIN+l/C 

[43]  XMIN+TZU-&C1 

[44]  r  YMINi  THE  MINIMUM  VALnE  OF  THE  COST  FUNCTION 

[45]  fl  FOR  SPECIFIC  T 

[46]  fl  XMIN:  THE  CORRESPONDING  AGE  REPLACEMENTTIME  (T) 

[47]  XA+XA,XMIN 

[48]  Y  A<rYA,YMIN 

[49]  fl  XAx  THE  VECTOR  OF  THE  AGE  REPLACEMENT  TIMES  (300) 

[50]  fl  YA:  THE  VECTOR  OF  THE  YMIN  (300) 

[51]  +(J<200)/LOOP2 

[52]  AXST+(+/XA)ipXA 

[53]  ACST+(+/YA)*pYA 

[54]  a  AXST:  THE  AVERAGE  VALUE  OF  THE  AGE  REPLACEMENT 

[55]  fl  TIMES  AFTER  300  REPETITIONS. 

[56]  fl  ACST:  THE  AVERAGE  VALUE  OF  THE  YMIN 

[57]  fl  AFTER  300  REPETITIONS. 

V 
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APPENDIX  C.  PROGRAM  AVERVVEIB 


c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  AVERWEIB 


t-,v*  ****** 


***  PURPOSE  :  This  program  calculates  the  average  costs  and  *** 
***  corresponding  average  replacement  times  by  using  the  *** 
***  sequential  estimating  procedure  for  the  weibull  distri-  *** 
***  bution.  The  program  also  calculates  the  mean  square  *** 
***  error  values  for  the  average  costs  and  the  average  rep-  *** 
*;f*  lacement  times  at  the  each  run.  *** 

**V r*  **  *  *  *  ***  *  *  ****  ******  ****  ****  *  *******  ***************  ********  * 


PARAMETER  (N=1000,  NREP=1000) 

INTEGER  I,  J,  K,  L,  M,  SEED,  DEL(N) 

REAL*8  X(N) ,  LAMHAT(N) ,  Z(N),  FISTAR(N),  XSTAR ,  SUMZ(N),  NUM(N) 
REAL*8  SUM(N) ,  SUM2(N),  AVG(N) ,  FSTAR,  MSEAGE(N) ,  Cl,  C2,  COSTST 
REAL*8  DEN(N) ,  ACOST(N),  SUM3(N),  SUM4(N),  AVACOS(N),  MSECOS(N) 


c 

*** 

The  following  real 

numbers  are  not  double  precision 

*** 

c 

*** 

because  the  random 

number  generator  subroutine  usee 

*** 

c 

*** 

single  precision 

only. 

*** 

REAL  U , LAM , AL 

c 

***; 

.-****************> 

fin 

m  mS+  *3+  l ^  ^  kt  ^  ^ 

?**** 

c 

*** 

KEY  VARIABLES 

*** 

c 

I  SEED 

The  random  number  seed 

*** 

c 

*** 

X 

Generated  system  lifetimes 

*** 

c 

*** 

LAM 

The  scale  parameter 

*** 

c 

*** 

AL 

The  shape  parameter 

*** 

c 

*** 

Cl 

Unplanned  replacement  cost 

*** 

c 

*3* 

C2 

Planned  replacement  cost 

*** 

c 

*** 

XSTAR 

A  constant  which  minimizes  the  cost 

*** 

c 

*** 

function  when  the  scale  parameter 

*** 

c 

.’-AJ. 

equals  1 

*** 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FSTAR 


COSTST 


ACOST 
FI STAR 


LAMHAT 


MSEAGE 


AVACOS 


MSECOS 


The  actual  age  replacement  time  (from  *** 
the  APL  program  Sim)  *** 
The  actua'i  cost  (from  the  APL  program  *** 
Sim)  *** 
The  return  generated  random  number  *** 
The  number  of  the  run  *** 
The  number  of  the  repetition  *** 
Indicates  the  summation  of  the  delta  *** 
values  to  calculate  the  average  costs  *** 
The  numerator  of  the  cost  function  *** 
which  given  in  the  Equation  (2.6)  *** 
The  denominator  of  the  cost  function  *** 
which  given  in  the  Equation  (2.7)  *** 
The  average  cost  *** 
The  average  age  replacement  time  *** 
The  minimum  value  of  the  X(i)  (genera-*** 
ted  system  lifetimes)  or  the  age  rep-  *** 
lacement  time  *** 
The  estimated  scale  parameter  *** 
The  average  value  of  the  age  replace-  *** 
ment  times  after  NREP  repetition  *** 
Mean  Square  Error  for  the  average  of  *** 
the  age  replacement  times  *** 
The  average  value  of  the  costs  after  *** 
NREP  repetition  *** 
Mean  Square  Error  for  the  average  of  *** 
the  costs  *** 


**********************************************-******************■ 
***  INITIALIZATION  *** 

SEED=16807 
LAM=0.  4557111695 
AL=1.  4D0 


Cl=10.  0D0 


c 


c 

c 

c 


C2=l. ODO 
FSTAR=0.  9169OD0 

***  From  the  main  lemma  xstar=f star/ lambda 

XSTAR=FSTAR/LAM 

C0STST=4.  0489 19D0 


DO  2  1=1, N 

SUM( I)=0.  ODO 
AVG(I)=0.  ODO 
MSE( I)=0.  ODO 
SUM2( I)=0.  ODO 
SUM3(I)=0.  ODO 
SUM4(I)=0.  ODO 
CONTINUE 


**************************************************************** 


***  CALCULATIONS 


£  ,Wc 


DO  3  L=1,NREP 

c 

CALL  EXCMS  ('FILEDEF  11  DISK  ALU  OUTPUT  A’) 

CALL  RANNUM  ( 1, SEED, 0.  0, 1.  0,0, U) 

C 

X( 1)=(( -1*L0G( 1-U))/(LAM**AL))**( 1/AL) 

Z(1)=X(1) 

LAMHATC 1)=1/X( 1) 

FISTARC 1)=LAMHAT( 1)*XSTAR 
SUMZ( 1)=X( 1)**AL 

C  ***  First  delta  value  is  1  (Failure  occur  before  the  time  *** 

C  ***  and  at  the  first  observation  there  is  no  comparison  for  *** 

C  ***  the  minimums.  Z(l)  =  X  (1) 

DEL( 1)=1 

C  ***  Calculation  of  the  average  cost  for  the  first  observation*** 

ktt  "X#/*  l  \  —  n  i 
ou;j(  i  )—\j  i 

DEN( 1)=Z( 1) 

AC0ST( 1)=NUM( 1)/DEN( 1) 
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CALL  RANNUM  ( 1 ,SEED,0. 0 , 1. 0 ,0 ,U) 


X(J)=((-l*LOG(l-U))/(LAM**AL))**(l/AL) 

***  Comparisons  to  find  the  Z  and  delta  values 
IF(FISTAR(J-1)  .LT.  X(J))  THEN 

SUMZ(J)=SUMZ(J-1)+(FISTAR(J-1)**AL) 

DEL( J)=DEL( J-l) 

Z( J)=FISTAR( J-l) 

ELSE 

SUMZ( J)=SUMZ( J- 1 )+( X( J)**AL) 

DEL( J ) =DEL( J - 1 ) + 1 
Z(J)=X(J) 

END  IF 

LAMHAK  J)=( DEL( J)/SUMZ( J) )**( 1/AL) 

FI STAR( J ) =LAHHAT ( J ) *XSTAR 

NUM( J )=( C 1*DEL( J ) )+( C2* ( J-DEL( J))) 

DEN( J)=DEN( J-l)+Z( J) 

ACOST(  J  )  =N”uM(  J  )  /DEN(  J  ) 

CONTINUE 
DO  10  K=1,N 

SUM(K)=SUM(K)+FISTAR(K) 

SUM2 ( K ) =SUM2 ( K ) +( F I STAR ( K ) -FSTAR ) **2 
SUM3(K)=SUM3(K)+AC0ST(K) 

SUN4(K)=SUM4(K)+( ACOST(K) -COSTST)**2 
CONTINUE 
CONTINUE 
DO  20  1=1, N 

AVG( I )=SUM( I ) /NREP 
MSEAGE( I )=SUM2( I ) /NREP 
AVAC0S( I )=SUM3( I ) /NREP 
MSECOSC I  )=SUM4(  I )  /NRF.P 

VRITE( 1 1 ,*)  AVG( I ) , AVACOSC I ) ,MSEAGE( I ) , MSECOSC I ) 


20  CONTINUE 
STOP 
END 

SUBROUTINE  RANNUM(DISTN,  SEED,  RPARM1,  RPARM2,  IPARM,  X) 

C  ***  THIS  SUBROUTINE  PROVIDES  AN  INTERFACE  WITH  THE  LLRANDOMII*** 

C  ***  ROUTINES  PROVIDED  IN  THE  NONIMSL  LIBRARY.  THE  PARAMETER  *** 

C  ***  REQUIREMENTS  AND  CALLING  PROCEDURES  ARE  AS  FOLLOWS:  *** 

0  VrVrVr  Jcirk 

C  ***  DISTN  =  DISTRIBUTION  TYPE  YOU  WANT  TO  SELECT  *** 

C  ***  AN  INTEGER  BETWEEN  1  AND  7  *** 

C  ***  SEED  =  THE  RANDOM  NUMBER  SEED  YOU  WISH  TO  USE  *** 

C  ***  RPARMl,  RPARM2 ,  AND  IPARM  ARE  REAL  AND  INTEGER  PARAMETERS*** 

C  ***  PASSED  TO  THE  ROUTINE  WITH  MEANINGS  WHICH  VARY  WITH  THE  *** 

C  ***  TypE  OF  DISTRIBUTION  YOU  DESIRE  *** 

C  ***  x  =  THE  RETURNED  RANDONM  NUMBER,  IT  IS  ALWAYS  REAL  *** 

C  *** 

C  ***  DISTRIBUTION  NUMBERS  AND  THE  ASSOCIATED  PARM  DEFINITIONS  *** 
0  *** 

c  ***  1- -UNIFORM  ON  THE  INTERVAL  RPARMl  TO  RPARM2  *** 

C  ***  2 --NORMAL  WITH  MEAN  RPARMl  AND  VARIANCE  RPARM2  *** 

C  ***  3- -EXPONENTIAL  WITH  RATE  RPARMl 

C  ***  4--C0UCHY  WITH  A  =  RPARMl  AND  B  =  RPARM2  *** 

C  ***  5 --GAMMA  WITH  SrlAPE  RPARM2  AND  RATE  RPARMl  *** 

c  ***  6--P0ISS0N  WITH  RATE  RPARMl  *** 

C  ***  7 --GEOMETRIC  WITH  P  =  RPARMl  *** 

REAL  RPARMl,  RPARM2 ,  X 
INTEGER  DISTN,  SEED,  IPARM,  N 

Q  ***  *** 

REAL  TEMP,  VARIAT(l) 

IF  (DISTN.  LE.O. OR. DISTN.  GT.  8)  THEN 

WRITE(10,  *)  'ILLEGAL  CALL  TO  RANNUM,  BAD  DISTN' 
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STOP 
END  IF 

C  Wt  *** 

GOTO  (10,  20,  30,  40,  50,  60,  70),  DISTN 
C  ***  *** 

c  ***  GENERATE  A  UNIFORM  BETWEEN  RPARM1  AND  RPARM2  *** 

10  CONTINUE 

C  **rr  *** 


IF  (RPARM1  -  RPARM2.EQ. 0)  THEN 

WRITECIO,  *)  'ILLEGAL  EQUAL  RPARMS  IN  RANNUM' 

STOP 

ENDIF 

IF  (RPARM1.GT. RPARM2)  THEN 
TEMP  =  RPARM1 
RPARM1  =  RPARM2 
RPARM2  =  TEMP 
ENDIF 

CALL  LRND( SEED ,  VARIAT,  1,  1,  0) 

VARIAT(l)  =  RPARM1  +  (RPARM2  -  RPARM1)  *  VARIAT(l) 

GOTO  99 

0  icicle 

C  GENERATE  A  NORMAL  WITH  MEAN  RPARM1  AND  STDDEV  RPARM2  *** 

20  CALL  LNORM( SEED ,  VARIAT,  1,  1,  0) 

WRITE(* ,  *)  'NORMAL  (0,  1)  ',  VARIAT(l) 

VARIAT(l)  =  (VARIAT(l)  *  RPARM2)  +  RPARM1 
GOTO  99 

Q  VoWr 

C  ***  GENERATE  AN  EXPONENTIAL  WITH  RATE  (1/MEAN)  RPARM1  *** 

30  CONTINUE 

IF  (RPARM1.EQ. 0)  THEN 

WRITE(10,  *)  'ILLEGAL  ZERO  RATE  IN  RANNUM* 

STOP 

ENDIF 

CALL  LEXFN(SEED,  VARIAT,  1,  1,  0) 
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VARIAT(l)  =  VARIAT(l)  /  RPARM1 
GOTO  99 

Q  ***  *** 

C  ***GENERATE  A  COUCHY  WITH  A  =  RPARM1  AND  B  =  RPARM2  *** 

40  CONTINUE 

IF  (RPARM2.LE.  0)  THEN 

WRITE(10,  *)  'ILLEGAL  COUCHY  SPREAD  IN  RANNUM,  B  =  ' ,RPARM2 
STOP 
ENDIF 

CALL  LCCHY(SEED,  VAR I AT,  1,  1,  0) 

VARIAT(l)  =  (VARIAT(l)  *  RPARM2)  +  RPARM1 
GOTO  99 

50  CONTINUE 

IF  (RPARM1.LE. 0)  THEN 

WRITE(10,  *)  ’ILLEGAL  NONPOSITIVE  GAMMA  RATE  IN  RANNUM* 

STOP  • 

ENDIF 

IF  (RPARM2.  LE.O)  THEN 

WRITE(10,  *)  'ILLEGAL  SHAPE  PARAMETER  IN  RANNUM' 

STOP 

ENDIF 

CALL  LGAMA(SEED,  VAR I AT,  1,  1,  0,  RPARM2) 

VARIAT(l)  =  VARIAT(l)  *  (1.0  /  RPARM1) 

GOTO  99 

60  CONTINUE 

IF  (RPARM1.LE.0)  THEN 

WRITE(10,  *)  'ILLEGAL  POISSON  RATE  IN  RANNUM’ 

STOP 

ENDIF 

CALL  LPOIS(SEED,  VARIAT,  1,  1,  0,  RPARM1) 

GOTO  99 

70  CONTINUE 

IF  (RPARM1.LE.O)  THEN 

WRITE(10,  *)  'ILLEGAL  GECM  PROB  IN  RANNUM' 
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STOP 

ENDIF 

CALL  LGEOM( SEED ,  VARIAT,  1,  1,  0,  RPARM1) 

GOTO  99 

***  *** 
CONTINUE 
X  =  VARIAT(l) 

END 
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APPENDIX  D.  PROGRAM  AVERGAM 


c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  AVERGAM 

ictrtcic  ItieirMrfch  *  VrVfVr 


***  PURPOSE  :  This  program  calculates  the  average  costs  and  *** 
***  corresponding  average  replacement  times  by  using  the  *** 
***  sequential  estimating  procedure  for  the  gamma  distribu-  *** 
***  tion.  The  program  also  calculates  the  mean  square  error*** 
***  values  for  the  average  costs  and  the  average  replacement*** 
***  times  at  the  each  run.  *** 
**************************************************************** 


PARAMETER  (N=1000,  NREP=1000) 

INTEGER  I,  J,  K,  L,  M,  ISEED,  DEL(N) ,  DELTA(N) 

REAL*8  X(N) ,  TEHAT(N),  Z(N),  FISTAR(N) ,  XSTAR,  COSTST 
REAL*8  SUM(N) ,  SUM2(N),  AVG(N) ,  FSTAR,  MSEAGE(N),  Cl,  C2 
REAL* 8  ACOST(N) ,  SUM3(N),  SUM4(N),  AVACOS(N),  MSECOS(N) 
REAL*8  NUMCO(N) ,  DENCO(N),  Tl,  T,  NUM,  A,  B,  C,  D 


c 

*** 

The  following  Real 

numbers  are  not  Double  precision 

*** 

c 

***  the  number  generator  subroutine  uses  single  precision 

REAL  U , TETA , P , RATE 

ititit 

c 

■sWrVf 

KEY  VARIABLES 

VnWr 

c 

ISEED  : 

The  random  number  seed 

*** 

c 

Vf** 

X  : 

Generated  system  lifetimes 

*** 

c 

J  0 

TETA  : 

The  scale  parameter 

c 

*** 

P  : 

The  shape  parameter 

*** 

c 

RATE  : 

Reciprocal  of  the  scale  parameter  to 

Vc  irk 

c 

*** 

use  the  subroutine  Rannum 

Vr-sV* 

c 

*** 

Cl  : 

Unplanned  replacement  cost 

*** 

c 

*** 

C2  : 

Planned  replacement  cost 

*** 

c 

*** 

XSTAR  : 

A  constant  which  minimizes  the  cost 

c 

function  when  the  scale  parameter 

*** 

c 

*** 

equals  1 

*** 
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c 

icicic 

FSTAR 

:  The  actual  age  replacement  time  (from 

*** 

c 

*** 

the  APL  program  Sim) 

*** 

c 

icicic 

COSTST 

:  The  actual  cost  (from  the  APL  program 

icicic 

c 

icicic 

Sim) 

*** 

c 

ft  ft* 

U 

:  The  return  generated  random  number 

*** 

c 

icicic 

N 

:  The  number  of  the  run 

*** 

c 

icicic 

NREP 

:  The  number  of  the  repetition 

id cic 

c 

*** 

DELTA 

:  Indicates  0  or  1.  If  failure  occurs 

icidc 

c 

icicic 

before  the  time  age,  then  equals  0. 

*** 

c 

icicic 

Otherwise  equals  1 

icicic 

c 

icicic 

DEL 

:  Indicates  the  summation  of  the  delta 

icicic 

c 

*** 

values  to  calculate  the  average  costs 

icicic 

c 

oWn’r 

NUMCO 

:  The  numerator  of  the  cost  function 

c 

which  given  in  the  Equation  (2*6) 

icicle 

c 

VnWc 

DENCO 

:  The  denominator  of  the  cost  function 

icicic 

c 

*** 

which  given  in  the  Equation  (2.7) 

icicic 

c 

*VoV 

ACOST 

:  The  average  cost 

icicic 

c 

FI STAR 

:  The  average  age  replacement  time 

icicic 

c 

^WcjV 

Z 

:  The  minimum  value  of  the  X(i)  (genera 

mirfcie 

c 

ted  system  lifetimes)  or  the  age  rep¬ 

4.^ 

c 

icicic 

lacement  time 

icicic 

c 

icicic 

TEHAT 

:  The  estimated  scale  parameter 

idcic 

c 

*** 

T1 

:  The  converged  scale  parameter  value 

icicic 

c 

at  the  E  step  of  the  EM  algorithm 

j+j.j- 

c 

VoWr 

T 

:  The  value  of  the  scale  parameter  at 

icicic 

c 

*** 

the  previous  calculation 

icicic 

c 

NUM 

:  The  numerator  of  the  Equation  4.  14 

icicic 

c 

VrVfVc 

A,  B,  C,  D 

:  To  determine  the  conditional  expecta¬ 

icicic 

c 

*** 

tion  value  from  Table  8 

icicic 

c 

*** 

AVG 

:  The  average  value  of  the  age  replace¬ 

icicic 

c 

ment  times  after  NREP  repetition 

icicic 

c 

VrVc’Vr 

MSEAGE 

:  Mean  Square  Error  for  the  average  of 

icicic 

c 

>V>Wf 

the  age  replacement  times 

icicic 

c 

*‘V:V 

AVACOS 

:  The  average  value  of  the  costs  after 
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c 

c 

c 

c 

c 

c 


c 


MSECOS 


NREP  repetition  *** 

:  Mean  Square  Error  for  the  average  of  *** 
***  the  costs  *** 

jVjWn'o'oV-iWf^VA^ViViVjVuViViWfjWdWfsWfiWciV*********************************** 

***  TMTTT  AT.T7ATTOM  * 


ISEED=16807 
TETA=2. 0/5.  0 
RATE=1. 0/TETA 
P=5.  0D0 
Cl=2. 0D0 
C2=l. 0D0 
FSTAR=1.  89690D0 

***  From  the  main  lemma  xstar=fstar*theta 


XSTAR=FSTAR*TETA 
C0STST=0.  933742D0 
DO  2  1=1, N 

SUM(I)=0.  0D0 
AVG( I )=0.  0D0 
MSE(I)=0. 0D0 
SUM2( I)=0.  0D0 
SUM3( I)=0.  0D0 
SUM4( I)=0.  0D0 
2  CONTINUE 

Q  it  V?  #*»«?*%  V*  Vc  Vc  *  ■/:  V:  Vr *  V?  MsMtitiei:' irtcirfrts&irtt’icic  itlrk' ititirtvh'liitiriiisiiiticisii'ii'kisiiitit' ititit 

C  ***  CALCULATIONS  *** 

Q  Vr  Vr  iV  ic  Vc  Vf  Vr  Vr  it  tV  tV  *  V  Vr  :V  Jc  V;  Vr  Vr  tV  V?  Vf  *  Vf  Vr  tV  Vr  Vr  V?  tV  Vr  rr  Yr  Yf  Yr  Yr  Vr  *  Yc  Yr  Yr  Yr  Yr  Y:  Yr  Yr  Yr  Yc  Yr  Yr  Yr  *  Yr  Yr  Yr  Yr  Yr  Yr  Yr  Yr  Yr  Yr  Yr  Yr  Yr 


DO  3  M=1 ,NREP 
C 

CALL  EXCMS  ( ’FILEDEF  15  DISK  P5C2  OUTPUT  A’) 
CALL  RANNUM  (5 , ISEED,RATE,P,0,U) 

C 

X(1)=DBLE(U) 
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Z(1)=X(1) 

C  ***  First  theta  value  from  the  Maximum  Likelihood  Estimation  *** 
TEHAT(1)=Z(1)/P 

C  ***  Use  the  main  lemma  to  find  the  estimated  age  replac.  time*** 

FISTAR( 1)=XSTAR/TEHAT( 1) 

C  ***  First  DELTA  value  is  1  (Failure  occur  before  time  age)  *** 

C  ***  and  at  the  first  observation  there  is  no  comparison  *** 

C  ***  Z(1)=X(1)  *** 

DEL(1)=1 
DELTA( 1)=1 

C  ***  Calculation  of  the  average  cost  for  the  first  observation*** 
NUMC0( 1)=C1 
DENC0( 1)=Z( 1) 

ACOST( 1)=NUMC0( 1)/DENC0( 1) 

DO  5  J=2 ,N 
C 

CALL  RANNUM  (5,ISEED,RATE,P,0,U) 

C 

X(J)=DBLE(U) 

C  ***  Comparison  to  find  the  Z  and  DELTA  values  *** 

IF(FISTAR( J-l)  .LT.  X(J))  THEN 
Z( J)=FISTAR( J-l) 

DEL( J)=DEL( J-l) 

DELTA( J)=0 
ELSE 

Z(J)=X(J) 

DEL( J)=DEL( J-l)+l 
DELTA( J)=l 
END  IF 
T=TEHAT( J-l) 

90  CONTINUE 

NUM=Z( 1) 

DO  100  L=2 , J 
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IF(DELTA(L)  . EQ.  1)  THEN 
C  ***  Observation  is  uncensored 

NUM=NUM+Z(L) 

ELSE 

C  ***  Observation  is  censored 

C  ***  Calculate  the  conditional  Expectation  E(X|X>Z)  values  *** 
C  ***  by  using  Table  8 

A=(Z(L)*Z( L) )+( 2*Z( L)*T)+( 2*T*T) 

B=(Z(L)*Z(L)*Z(L) )+( 3*T*A) 
C=(Z(L)*Z(L)*Z(L)*Z(L))+(4*T*B) 

D=( Z(L)*Z(L)*Z(L)*Z(L) )+(5*T*C) 

NUM=NUM+(D/C) 

END  IF 

100  CONTINUE 

T1  =  NUM  /  (J*P) 

IF  (ABS(Tl-T)  .LT.  0.001)  GO  TO  200 

C  ***  This  if  statement  is  to  check  stopping  criteria  *** 

C  ***  If  this  satisfied  we  can  accept  T1  converged  to  the  theta*** 

T=T1 

GO  TO  90 
200  CONTINUE 

TEHAT( J)=T1 

FISTAR( J)=XSTAR/TEHAT( J) 

NUMC0( J )=( C 1*DEL( J) )+( C2*( J-DEL( J) ) ) 

DENC0( J)=DENC0( J- 1 )+Z( J) 

AC0ST( J)=NUMC0( J)/DENC0( J) 

5  CONTINUE 

DO  10  K=1 ,N 

SUM(K)=SUM(K)+FISTAR(K) 

SUM2(K)=SUM2(K)+(FISTAR(K) -FSTAR)**2 
SUM3(K)=SUM3(K)+AC0ST(K) 

SUM4(K)=SUM4(K)+( ACOST(K) -C0STST)**2 
10  CONTINUE 

3  CONTINUE 
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DO  30  1=1, N 

AVG( I )=SUM( I ) /NREP 
MSEAGE( I )=SUM2( I)/NREP 
AVACOS( I )=SUM3( I ) /NREP 
MSECOS( I)=SUM4( IJ/NREP 

WRITE ( 15 , 300 )  AVG( I ) , AVACOS( I ) ,MSEAGE( I ) , MSECOS( I ) 
300  FORMAT  (F10.  7,'  '.F10.7,'  ',F10.7,' 

30  CONTINUE 
STOP 
END 


'  ,F10.  7) 
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